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Abstract 

We study wave propagation in networks of coupled cells which can be- 
have as excitable or self-oscillatory media. For excitable media, an asymp- 
totic construction of wave trains is presented. This construction predicts 
their shape and speed, as well as the critical coupling and the critical sepa- 
ration of time scales for propagation failure. It describes stable wave train 
generation by repeated firing at a boundary. In self-oscillatory media, 
wave trains persist but synchronization phenomena arise. An equation 
describing the evolution of the oscillator phases is derived. 

1 Introduction 

Understanding wave propagation and self-oscillations in discrete media is impor- 
tant to predict the behavior of many biological and physical systems composed of 
smaller interacting units, such as cells, nodes, atoms, wells... Examples abound: 
propagation of impulses through myelinated axons [2SI or cardiac tissue [2], 
voltage oscillations in muscle fibers |22j. charge transport in semiconductor su- 
perlattices crystal growth and interface motion in crystalline materials [H], 
etc. 

In this paper, we consider an infinite chain of diffusively coupled cells de- 
scribed by a fast variable v n and a slow variable w n : 

dv 

— = D(v n+1 - 2v n + + f(v n ,w n ), (1) 



dt 



= Xg(v n ,w n ). (2) 



Typically, v n (t) represents the ratio of the membrane potential to some reference 
equilibrium potential and w n (t) the fraction of open channels for some type of 
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ions at the n-th site. The parameter D > measures the strength of the 
coupling and A > represents the ratio of the time scales for v n and w n . A 
must be small enough for the two variables to evolve in different time scales. 
We assume that the nullclines g(v,w) = = f(v,w) intersect at one point, the 
equilibrium of the system. For w fixed, the source f(v,w) is bistable. System 
CJ-® displays excitable or self-oscillatory behavior depending on whether the 
equilibrium state is stable or unstable. 

In the excitable case, periodic disturbances in the medium lead to propaga- 
tion of signals in the form of wave trains. To obtain analytic information on 
those waves, it is quite tempting to replace the discrete system Q-© by its 
continuum limit D — ► oo: 

dv d 2 v . . dw . . . . 

oi = W + fiv > w)t -m =X9 ^ (3) 

There is a good mathematical reason not to do that: propagation failure. The 
discrete system has a threshold coupling for wave propagation whereas the con- 
tinuous model allows wave propagation at all coupling strengths. Constructing 
wave train solutions of discrete systems is not easy. An asymptotic understand- 
ing of these solutions and their numerical construction have to go hand in hand. 
In an infinite chain, most initial conditions will either fail to propagate or evolve 
into a solitary pulse. Only very particular initial conditions evolve into a wave 
train, as we will show later in this paper after explaining the asymptotic con- 
struction thereof. Simpler pulse solutions were constructed asymptotically in [5] 
for the discrete FitzHugh-Nagumo model. For continuous models, asymptotic 
descriptions of pulse waves are given in Unlike continuous model equations, 
travelling waves in discrete media do not solve a system of ordinary differential 
equations. Instead, a more complex system of differential-difference equations 
needs to be analyzed. Rigorous results are only available for particular sources 
and traveling pulses. Tonnelier found explicit expressions of pulses for piece- 
wise linear sources. Existence of discrete pulses has been proved for a restricted 
class of nonlinear sources |1 lj . Conditions for propagation failure were discussed 
in Q]03|Ii]. 

In this paper, we give an asymptotic construction of wave trains in discrete 
excitable systems. Solitary pulses are a particular case of a wave train with 
infinite spatial period and maximum speed. We study the dynamic equations 
f° r general smooth nonlinear source terms. Experimental values of D 
(for instance in nerve models) are typically small. Thus, we focus in the highly 
discrete case: D <C 1. We find a one-parameter family of stable wave trains 
and we predict their speed, period and shape, as well as the critical coupling 
for propagation failure. Near the critical coupling, the profiles develop steps as 
shown by Figure E£b) and Figure 03 Propagation is 'saltatory': it proceeds by 
a sequence of successive jumps. The mechanism for propagation failure as the 
coupling decreases is reminiscent of the depinning transitions for wave fronts 
described in Propagation also fails when the parameter A characterizing 
the ratio of time scales is no longer small. Our asymptotic construction yields 
sharp upper an lower bounds on the critical value of A below which propaga- 
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tion fails, X C (D), as a function of the coupling D. For D small, we also obtain 
an asymptotic prediction which describes accurately numerical measurements. 
Technically, we have greatly improved the rough bounds contained in de- 
rived under very restrictive assumptions on the nonlinearities. Bounds on the 
critical separation of scales are unknown in the continuum case. 

When the medium is self-oscillatory, the mechanism leading to periodic wave 
train generation is completely different. Spatially uniform profiles corresponding 
to stable limit cycles of the system: 

-=f(v,w),^ = Hv,w), (4) 

play a key role and introduce new features such as synchronization. For small 
couplings, we derive an equation for the time evolution of oscillator phases. 

The paper is organized as follows. In section |3 we give an asymptotic de- 
scription of periodic wave trains for excitable systems. Section [3] discusses their 
numerical stability and selection via initial and boundary conditions. We use 
the FitzHugh-Nagumo (FHN) and the more complicated Morris-Lecar (ML) 
dynamics as examples in our numerical calculations. Section 0] analyzes the im- 
pact of discreteness on propagation failure. We predict the critical coupling and 
the critical separation of time scales for propagation failure as a function of the 
coupling. Sectional is devoted to self-oscillatory behavior and synchronization. 
Section El contains our conclusions. 



2 Wave trains and pulses in excitable media 

Figure n depicts the nullclines for system Q-© and the sources we have used 
in the numerical tests. We assume that the source f(v, w) is a cubic function of 
v with three zeros z\{w) < z%[w) < 23(11)), the first and the third of which are 
stable for the reduced equation Q with w n — w. The nullclines must intersect 
in the first or third branch of the cubic to guarantee that ©-(0) describe an 
excitable system. Intersections of the nullclines in the middle branch of the 
cubic correspond to unstable stationary states and give rise to self-oscillatory 
dynamics. To fix ideas, we assume in this section that the nullclines intersect 
at a point (v Q ,w°) located in the first branch of the cubic. 

Experimental values of D (for instance in nerve models) are typically small. 
Thus, we focus in highly discrete problems with D « 1. As shown in Figure 
12 a wave train resembles a sequence of pulses. After a short transient the 
system relaxes to a traveling wave: v n (t) — V(n — ct), w n (t) = W(n — ct), in 
which V and W are periodic functions with the same period L. This defines the 
wave profiles V{z) = v n { 11 ^-) and W(z) = w„(^^), z = n — ct. Neighboring 
points describe the same trajectory with a time delay -: v n +i(t) = v n (t — -), 
w n+ i(t) = w n (t - i). 

The profiles and speeds solve an eigenvalue problem for a system of diffcrcntial- 
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Figure 1: (a) FHN with a = 0.3 and b = 0.5 (solid), b = 10 (dotted); (b) ML 
with I = 0.0625 (dashed), 0.375 (dotted). 



difference equations: 

-cV z (z)=D(V(z+l)-2V(z) + V(z-l)) + f(V(z),W(z)), 

-cW z {z) = Xg(V(z),W(z)), [0) 

with V(z + L) = V(z), W(z + L) = W(z). Equations J5J are hard to solve 
except in very particular situations, such as piecewise linear / and linear g 
|26| . Is there an alternative way to extract practical information on the waves? 
We show below how to use asymptotic methods to predict the ranges of the 
parameters for which we expect stable wave trains to be formed, the shape of 
their profiles, their spatial period L and speed c. 

2.1 The reduced bistable equation 

As FigureElshows, a wave train for v n is a sequence of sharp interfaces separated 
by smooth transitions. At the interfaces, the recovery w n is almost constant. 
Thus, the first step in our asymptotic construction is the study of wave front 
propagation in the bistable equation: 

dv 

= D(v n+1 - 2v n + v n _i) + f(v n , w). (6) 

Choosing w in a neighborhood (u>_, w + ) of the equilibrium value w°, the source 
f(v, w) has three zeros z\{w) < Z2(w) < z^(w), the first and the third of which 
are stable. The parameter w controls the degree of asymmetry of the source, as 
shown in Figure |3Jb) . In particular, it changes the ratio of the area A12 above 
/ between z\ and zi and the area A23 below / between Z2 and Z3. When the 
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cubic source f(v, w) is asymmetric enough, depending on D, wave fronts v n — 
V(n— ct) joining the two stable zeros Z\ O) and z 3 {w) propagate [71151 1111011 127) . 
Generically, there are two threshold values w c - (D) ,w c + (D) such that: 

• If w € (w- ,w c -(D)), wave fronts decreasing from z^(w) to z\{w) (resp. 
increasing from zi(w) to z^fw)) move to the right (resp. left) with speed 
c = c+(w, D). 

• If w S (w c+ (D),w + ), wave fronts increasing from Zi(w) to Zs(w) (resp. 
decreasing from z^(w) to move to the right (resp. left) with speed 
c = C-(w, D). 

In the first case, A23 > A 12 . In the second, A 2 z < A\ 2 . If w € {w c -(D),w c+ (D)}, 
fronts are pinned and c = 0. The size of this pinning interval grows as D 
decreases. 

Let us recall some facts about the structure of wave front solutions, which 
will be needed in Section 0] For small D, the stable stationary front s n with 
w = w c -(D) has the following structure: two essentially constant tails s n ~ 
zz{w), n < 0, and s n ~ Zi(w), n > 0, joined by an active point sq, which takes 
values in the interval (zi(w), z 2 {w)). When w is near the threshold w c -(D), the 
travelling wave front profiles are staircase-like. While the wave fronts move at 
a constant speed c, the motion of points v n (t) belonging to these wave fronts 
is saltatory: their spatial profiles v n stay near the shifted static configurations 
s„_;, I £ Z, for a time -, until the l-th point Vi(t) (the active point) jumps from 
So to zz{w) and vi+±(t) (the next active point) jumps from Z\ (w) to sq, making 
the front advance one position. Increasing D, we find a similar situation with an 
increasing number of active points filling the gap between the constant tails in 
the spatial profiles and an increasing number of steps in the temporal profiles. 
As w grows, the profiles and their motion become smoother, see [51ll2|. 

The speed of the wave fronts can be predicted in two limits. In the depinning 
transition w — > w c ±(D), we find c w sign(u> — w c ±{D))^f a/3(w — w c ±(D))/ir, 
see 8 for details. For highly discrete problems with D« 1, fronts can move 
only when w ~ w c ±(D) and the approximation for the speed provided by the 
analysis of the depinning transition always holds. This scaling changes when / 
is a piecewise linear function, as shown in |12| . In the continuum limit fl>l, 
the wave front velocities are a correction of the wave front velocities cq for a 
spatially continuous problem c w y/Dco(l - see |T5J for details. Out of 

these two limits, characterized by either small or large speeds, we compute the 
speeds by solving JHJ numerically. 

The picture of wave front propagation we have described holds for any 
smooth source / when D is small enough. For larger D, we may encounter 
sources / for which the pinning interval shrinks to a point 0. Finding a com- 
plete characterization of the sources for which there is no pinning is an unsolved 
problem. These anomalous cases are excluded from our analysis. 
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Figure 2: Wave train for FHN: (a) Temporal profiles of v n (t) (solid) and w n (t) 
(dashed-dotted), compared to their asymptotic reconstruction (dotted); (b) Spa- 
tial profiles of v n (circles) and w n (triangles). Parameter values are a = O.f, 
b = 0.5, D = 0.02, A = 0.01, w = 0.002. 



2.2 Asymptotic construction of wave trains 

Having described the reduced dynamics of wave front solutions of bistable equa- 
tions in Section f2.il we are ready to predict asymptotically the speed and shape 
of wave trains for m the limit A — * + . 

Discrete wave trains have the form v n (t) = V(n — ct), w n (t) = W(n — ct). 

V and W are real junctions of a real variable z = n — ct. Note that z takes on 
continuous values and not discrete values like the variable n. Thus V and W are 
perfectly determined once their speed c, period L and continuous profiles V(z) 
and W(z), z G [0, L], are known. If we fix n, the trajectories v n (t) = V(n — ct), 
w n(t) = W(n — ct) are periodic in the continuous variable t and have a well 
defined temporal period L/c. Thus, once we compute the speed c and the 
temporal profiles v n (t) , w n (t) as functions oft for fixed n, the wave train profiles 

V and W and their period L are known. If we fix t and vary the discrete spatial 
variable n, the situation is different. The spatial structure of the wave train at 
each fixed time t is given by v n — V(n — ct) and w n — W(n — ct). The points 
contained in a period satisfy < n — ct < L, that is, ct < n < L + ct. As 
time grows, the discrete points 'travel' along the continuous wave profile and 
are transferred from one period to the next. The number of integers n one can 
fit in an interval of length L is the integer part of L. 

There are no rigorous existence results for wave train solutions of JU-0- 
Assuming that wave trains do exist, we shall use the separation of time scales 
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A to find an approximate reconstruction of the temporal wave train profiles by 
matched asymptotic expansions at zeroth order as A — * + . A description of 
matched asymptotic expansions can be found in [18] or in the more elementary 
reference [3]. 

For D small, numerical solutions show that wave train trajectories should 
be composed of regions of smooth variation of (v n , w n ) on the slow time scale 
T = Xt, separated by sharp interfaces in which v n varies rapidly in the fast time 
scale t, see Figure [21 The reference time in our construction of trajectories will 
be the slow time T. 

For n fixed, we may distinguish four stages in one time period of the wave 
train solutions (v n (T; X),w n (T; A)). At each stage, a reduced description holds: 

• Leading edge (I) in Figure [21 (a) . The variable v n evolves abruptly in the 
fast time scale t = T ~ T ° £ (—00,00). To leading order: 

^ = D(v n+l -2v n + v n -i) + f(v n ,w n ), *6t=0. (7) 

Thus, w n = w is a constant and v n is the temporal profile of the wave 
front solution of JfjJ which joins Zx(w) and 23 (w). This front propagates 
with a speed c+(u>, D) when w £ (w-,w c -(D)), 

• Peak. This is region (II) in Figure [21 (a). Here v n varies smoothly on the 
slow time scale T. Thus, we may set D = 0, obtaining the reduced problem 
f(v n ,w n ) = 0, = g(v n ,w n ). Matching the peak to the leading edge 
at zeroth order, we see that v n must lie in the third branch of roots of the 
cubic function: v n — Z3(w n ). Then, for Tq < T < Ti(w,w): 

= g(z 3 (w n ),w n ) : (8) 

T 1 ( W ,uJ)-T + / f (9) 

in which w n (To) — w and w n (Ti) = w (to be calculated later). 

The matching conditions at T = To for the reduced descriptions of the 
peak and the leading front are 1 : 

zM H) (T)) - vlP(^) « 1, wi H \T) - w « 1, 

as A — * 0+ in the overlap region A<T-To<l. The superscripts (I), 
(II) ... refer to the region and the reduced description we are using. The 
uniform approximations of the profiles in stages (I)- (II) are: 

v umf = „W("i ) + Z3 {w (iD ( T ))-z 3 (w), =wk U) (T). 



1 /(e) <C as e — > is equivalent to / = o(g) and it means lim e ^o//9 = 0, see [3|. 
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Trailing wave front. The peak is followed by another region of fast vari- 
ation of v n , corresponding to (III) in Figure [5] (a). Again, Q holds to 
leading order, on the new time scale t = T ~ x Tl S (—00, 00). Thus, w n = w 
is a constant and v n is the temporal profile of the wave front solution of 
lj?Jl with w n — w which joins z 3 (w) and Z\(w). This front propagates with 
a speed C-(w,D) when w S (w c+ (D),w+). If the wave train is to move 
rigidly, w must be selected in such a way that the trailing front travels 
with speed c_(W, D) = c + (w, D) = c. 

The matching conditions for the reduced descriptions of the peak and the 
trailing front at T = T\ are: 

z 3 (wi H \T)) - vk IU) (^) « 1, -w n II] (T) + w « 1, 

as A — > 0+ in the overlap region A <C T± — T <C 1. The uniform approxi- 
mations of the profiles in regions (I)- (III) are: 

v umf = W (^)(Z_Z1) + + z,{w^\T)) - z 3 (w) - z 3 (w), 

= wi /7) (T). 

TaiZ. The last part of the period is region (IV) in Figure |5] (a). Here, v n 
varies smoothly on the slow time scale T, and D can be neglected. The 
reduced problem is f{v n ,w n ) = 0, ^f- — g(v n ,w n ). Matching the tail to 
the trailing wave front at zeroth order, we see that v n must lie in the first 
branch of roots: v n = zi(w n ). Then, for T% <T < T2(w,w): 

= 9(zi(w n ), w n ), (10) 
ds 



T 2 {w,w)=T 1+ / - — — — - , (11) 

Jw g(zi(s),s) 

in which w n (Ti) = W and w n {T2) = w = w n (To), if the wave train is to 
move rigidly. In this way, we conclude the description of one temporal pe- 
riod and guarantee that the end of its tail can be matched to the following 
(identical) period. Notice that the integral ljll|l diverges when w > u>° 
due to a singularity at the equilibrium value w . Thus, the range of w for 
which we can find wave trains is restricted to (w° ,w c -(D)). 

The matching conditions for the reduced descriptions of the tail and the 
trailing front at T — T\ are: 

Vn HI \^) ~ ZM IV \T)) « 1, -W n IV \T) + W « 1, 

as A — > 0+ in the overlap region A -C T — T\ <C 1. The uniform approxi- 
mations of the profiles in regions (III)- (IV) are: 

y umf = v ini) { T^ L) + Zl(w f)( T) ) _ w umf = ^(T). 

In a similar way, the end of the tail should be matched to the leading edge 
(region (V)) of the next period at T^. 
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A uniform approximation to the temporal profile in a whole period is ob- 
tained as follows. Let x( a i be the indicator function, equal to 1 if a < T < b 
and otherwise. Then, for T <T <T 2 : 

v um f{T) = [ w (^) ( z^ )+w a) ( z^ )+Z 3 (w (^)(r))-z3H-z 3 (w)]x(T ,T 1 ) 

+ [^ 7) (^) + ^ ) (^) + z 1 ( w r ) (T))-, 1 H-z 1 (Tiy)] X (T 1 ,T 2 ), 
w um 1{T) = W W {T ) X (T , Tx) + w^iTMT,, T 2 ). 

This construction yields a family of wave trains parametrized by the values 
w G (w , w c -(D)) for which there exists w = w(w) € (w c +(D),w+) satisfying 
c_ (mJ, D) = c + (w, D) = c. 

In Figure |2] (a) , we have superimposed the approximate temporal profile 
provided by our matched asymptotic expansions and the numerically calculated 
profile. They are indistinguishable. Notice that, for small D, the leading and 
trailing wave fronts can be accurately computed using one active point (see 
section 2.1 or [8]): 

D(z 3 (w) - 2v„ + zi(w)) + f(v n ,w), (12) 

D(z 3 (w) - 2v n + Zl (w)) + f(v n ,w). (13) 

In the true time scale t = ?j the temporal period of the wave train is asymptotic 
to (T2(w) — To) /A (the duration of the fast stages (I) and (III) is ignored). 

Having calculated the temporal profiles v n (t),w n (t) of a wave train and its 
velocity c, the wave train profiles V, W are known: V{z) — v n ( 11 ^) and W(z) — 
w„( Ii ^ £ ), z = n — ct. Their period is now asymptotic to c + (w)(T2(w) — T )/A. 
The number of points L(w) in the wave length of the wave train is essentially 
the integer part of the wave length. Figurc[5](b) illustrates that, at a fixed time 
t: 

• the leading front of one spatial period is located at a point no, where 

W n[) (t) ~ W, 

• its trailing front is located at a point m, such that w ni (t) ~ w, 

• the end of our spatial period and the beginning of the next one is located 
at a point n 2 , where w n2 (t) ~ in. 

Then, L(w) = no — n 2 is approximately the integer part of c+ (w) (T 2 (w) — To) /A. 
Notice that the number of points inside the leading and trailing fronts of each 
period can be ignored in the discrete limit D<1. The number of points in the 
peak Li(w) = no — n\ is asymptotic to the integer part of c+(u>) \T\{w) — To) /A 
and the number of points in the tail L 2 (w) = n\ — n 2 is asymptotic to the integer 
part of c+(w)(T 2 (w) — Ti(w))/X. Our construction is consistent only when 
L\(w) > 1. This provides an estimate of the critical value A C (D) above which 
propagation fails. Nevertheless, near the thresholds for propagation failure, the 
approximation described in this section has to be corrected. This is explained 
in more detail in Section 0] 



dv n 

dt 

dv n 

dt 
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2.3 Solitary pulses 

A wave train becomes a traveling pulse when its spatial period tends to infinity, 
which occurs as w — > w°. Pulses are the fastest waves of the family. In the 
asymptotic description of the temporal profile of a pulse we distinguish five 
regions: 

• In front of the pulse, the profile is at equilibrium: v n ~ v , w n ~ w . 

• The leading edge of the pulse is a wave front solution of (JJJ with w n — w° , 
which joins Z\(w ) = v° and 23(111°). This front propagates with a definite 
speed c = c + (w° , D) . 

• In the transition between fronts, v n — Z3(w n ), = g{zz{w n ),w n ) and 
w n evolves from w° at T to w° at T\ . 

• The trailing wave front is the wave front solution of (JZJ) with w n = w°, 
which joins Z3(w°) and Zi(Hr). w° is now selected in such a way that this 
front travels with speed c = c_(TtJ°, D) — c+(w°, D). 

• In the pulse tail, v n = zi(w n ), ^§f- = g(zi(w n ), w n ) and w n evolves from 
w° to w°. 

The number of points between fronts Li(w°) is again asymptotic to the 
integer part of of c + (w°)(Ti(w°) — T )/X. The number of points in the tail is 

infinite, since the integral JLq g ( z f* s ) s ) diverges. However, we can predict how 
long does it take for the tail to get sufficiently close to v° — zi(w°) by using 
f—o 6 , d , s , r , with e > instead. 



3 Periodic firing of waves at boundaries 

When an excitable medium is periodically excited, we expect propagation of 
signals in form of wave trains. We will see here that the periodic wave train 
solutions constructed in Section fallow to describe propagation due to periodic 
firing at a boundary in the spatially discrete FHN and ML models. 



3.1 FitzHugh-Nagumo 

The FitzHugh-Nagumo (FHN) model is a crude simplification of Q-© that 
sets f(y,w) = h(v) — w with g linear. Typically, h is chosen to be a cubic 
polynomial, say h(v) = v(v — a) (2 — v). We get: 

dv 

-£=D(v n+1 -2v n +v n _ 1 )+v n (v n -a)(2-v n )-w n , (14) 

—7^- = A(«„ - bw n ). (15) 
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Figure 3: For FHN: (a) Critical curves for cubic shape; (b) Change ol symmetry 
in the source as w increases from to u>+; (c) Speed curves when D = 0.1 for 
a = 0.55 (solid), a = 0.3 (dashed), a = 0.1 (dashed-dotted); (d) Dispersion law 
L(c) for D = 0.1 and A = 0.05. 



When a £ (0, 1), this system is excitable provided b is not too large, as illus- 
trated in Figure^a). The unique stationary state is (v°, w°) = (0, 0). As Figure 
|3Jb) shows, the source v(v — a) (2 — v) — wis bistable when w £ (iw_(a), w+(a)). 
The curves w-(a), w+(a) are plotted in Figure |SJa). The reduced equation l|14fl 
with w n = w has decreasing (resp. increasing) traveling wave front solutions 
moving with speed c+(w) (resp. C-(w)) when w e (w-(a), w c -(a, D)) (resp. 
w e (w c +(a, D), w+(a))). FigureOIc) depicts the speed curves c(w) for different 
values of a. The left branches represent c+(w) and the right branches c_(u>). 
The horizontal dotted lines join c+(0) with c_(w°), selecting the value uf for 
the trailing front in the construction of pulses. 

Wave trains are found when w £ (0, w c —(a, D)). They are slower than 
pulses. Figure0{a) shows a wave train for a = 0.3, w = 0.02 and A = 0.01. It 
has been generated by solving (|14[l - i|15|) with initial data v n (0) = 2, n = 1, 10, 
fn(0) = 0, n > 10, uvt(0) = 0.01 and a periodic excitation at the left boundary. 
More precisely, we have used boundary conditions of the form: 

( ] _ f 2 t£ [0,pT] + kT, 

with k G JV, < p < 1. Our theory predicts c = 0.18, L\ = 9 and L 2 = 45, 
close to the numerically measured values c = 0.175, L\ = 10, L2 = 46. For 
w = 0.002, our theory predicts c = 0.2. When A = 0.05, Li = 2 and L 2 = 14. 
Numerically, we find c = 0.16, Li = 1 - 2, £ 2 = 13 - 14. 
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Figure 4: (a) Wave train for FHN when D = 0.1, a = 0.3 A = 10~ 2 , w = 0.02, 
c = 0.18; (b) Wave train for ML with D = 1, A = 10~ 3 , w = 0.045, c = 0.54. 



We have plotted the dispersion relation L(c) (L(w) against c+(w)) for our 
asymptotic family of wave trains in Figure Old). Note that L'(c) > 0. Our nu- 
merical tests with periodic boundary conditions starting from initial conditions 
which are near the shape asymptotically predicted for these wave trains suggest 
their stability. For continuous systems like (J3J), Maginu ^5] proved stability of 
wave trains under the condition L'(c) > 0. Our numerical experiments suggest 
that a similar stability criterion could be established for wave trains in discrete 
systems. 



3.2 Morris-Lecar dynamics 

The non dimensional Morris-Lecar (ML) model [22] is: 

dv 

= D{v n +i - 2v n + «„_!) + f(v n ,w n ) - 21, (17) 

— = A cosh( — ) [1 + tanh( — ) - 2w n \ , (18) 

where the index n denotes the n-th site and: 

f(v, w) = 2w(v - V K ) + 2g L (v - V L ) + g Ca [l + tanh(^— -)] (v - 1). (19) 

v n is the ratio of membrane potential to a reference potential and w n is the 
fraction of open K + channels. The time scale is , ~g~K being the K + conduc- 
tance and C m the membrane capacitance. System (Il7ll - <|18|) is a reduced version 
of the full Morris-Lecar model which involves one more fast variable m n . 
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Figure 5: For ML: (a) Change of shape in f(v,w); (b) Pinning thresholds 
w c -(D) and w c +(D); (c) Speed curves C-(w) and c + (w) with D = 1; (d) 
Dispersion law L(c) for D = 1 and A = 10~ 3 . 



Typical values for the parameters can be taken from experiments |22|. In our 
numerical tests we have used: 
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System (I17|) - (I18|I exhibits a rich dynamical behavior depending on its null- 
clines. Figure ^b) shows two possibilities. For / = 0.375, the system develops 
self-oscillations. This behavior will be discussed in Section [S] For / = 0.0625, 
the system is excitable. Figure OJa) illustrates the change of symmetry in the 
source f(v,w) as w increases from w- to w+. The thresholds for wave front 
propagation w c -{D) and w c+ (D) are depicted in Figure GJb). Note that the 
size of the pinning interval grows as D decreases. The curves C-(w) and c + (w) 
representing the speeds of the wave front solutions of the reduced bistable equa- 
tion for this problem with D fixed are shown in Figure 03c)- We have plotted 
the dispersion relation L(c) (L(w) against c + (w)) for our asymptotic family of 
wave trains in Figure EJd). Again, L'(c) > and this branch of wave trains is 
expected to be stable for the dynamics lfT7|) - ((TH|l . 

We have generated a variety of wave trains using the asymptotic prediction 
as initial datum and periodic boundary conditions. These numerical solutions 
agree reasonably well with the asymptotic description given in Sections 12.21 
Figure 0Jb) shows a wave train for D = 1, A = 10~ 3 , w = 0.045, c = 0.54. 
For w — 0.045, the predicted speed is c = 0.54. When A = 10 -3 , we predict 
a peak width L\ ~ 19 and a tail width L2 ~ 256. The numerically measured 
parameters are c ~ 0.55, L\ ~ 18 and L2 ~ 250. 
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The basin of attraction of wave trains is rather small and their existence 
depends crucially on the boundary conditions we employ There are no wave 
trains if we replace the periodic boundary conditions by Dirichlet or Neumann 
boundary conditions. With periodic boundary conditions, wave trains fail to 
be formed if w n (0) is near the equilibrium value w°. In such a case, the initial 
condition evolves toward a solitary pulse instead of a wave train. 

Pulses are more robust. They can be generated using different initial data 
and boundary conditions. For instance, we can chose step-like initial data v n (0), 
climbing from v° to a value close to 23(117°). For w n (0) we may take a pertur- 
bation of w° , as long as w n (0) is near the value w° in the trailing front of the 
step. These initial conditions evolve into pulses, regardless of whether we are 
using Dirichlet or periodic boundary conditions. When D — 0.1, our asymp- 
totic construction predicts pulses with speed c ~ 0.075. Setting A = 10~ 4 , the 
predicted peak width is L\ ~ 35. The tail reaches the value w° — e, e = 10~ 2 , for 
L e ~ 455. Again, these values describe accurately the numerically constructed 
pulses. 

We observe that the profiles of wave trains and pulses differ slightly de- 
pending on whether their speeds are large or small. Fast waves have smooth 
temporal profiles and advance smoothly. Slow waves develop a sequence of steps 
in the leading and trailing fronts. Their motion is 'saltatory'. Figure HH shows 
the spatial and temporal profiles of a slow wave train when D = 1, A = 10~ 5 , 
w = 0.108 and c = 0.04. The values w and w are close to the thresholds u> c _ (D) 
and w c +(D), that is, near the depinning transitions |S] for the bistable equation 
(0. This is the reason for the appearance of the steps that can be appreciated 
in Figure HJb). Wave trains and pulses are always 'slow' when D is near the 
threshold D c for propagation failure. For our parameter values, w° is near the 
threshold w c _(D) when D ~ 0.076, and the profile of the pulses and wave trains 
are staircase like. 

4 Propagation failure 

In both discrete and continuous models, wave trains and pulses disappear as 
A increases beyond certain thresholds. However, the way in which these waves 
disappear can be rather different. In continuous models, the height of the peak 
of a pulse or a wave train decreases as A increases and, eventually, waves cease to 
propagate. In discrete models, new modes of propagation failure arise depend- 
ing on the value of D. Unlike the continuous case, waves also fail to propagate 
when D is smaller than a critical value D c . Figure [7| depicts the critical sepa- 
ration of time scales A C (D) for FHN as a function of D. We distinguish three 
regions: small D (highly discrete limit), large D (continuum limit) and inter- 
mediate D. For the parameter values selected in Figures [71101 these regions are 
approximately D c < D < 0.01, 0.01 < D < 0.3 and 0.3 < D. FiguresEI(b) and 
(c) zoom in the highly discrete and the intermediate regions. 

Figure |H1 illustrates the shape of pulses for FHN when D is near D c . For A 
close to \ C (D), the height of the pulses has not decreased and remains near its 
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Figure 6: Slow wave trains for Morris-Lecar: (a) spatial profile at a fixed time 
t, (b) detail of the steps in the temporal profile. Parameter values are D = 1, 
A = 1(T 5 , w = 0.108 and c = 0.04. 
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Figure 7: Numerically measured X C (D) (asterisks) compared to: (a) the upper 
bound X C (D) + (circles) and the lower bound A C (Z?)~ (triangles); (b) and (c), 
the upper bound \ C (D) + (circles), the numerically measured Ao(-D) (squares) 
and its asymptotic prediction (|25(l (triangles). Parameter values for FHN are 
a = 0.1, b = 0.5, k = 5. 
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Figure 8: Pulses near D c and A C (Z?): (a) Spatial profile and (b) Temporal 
profile (solid line) when D = 0.004 and A = 0.0054; (c) Spatial profile and (d) 
Temporal profile (solid line) when D = 0.009 and A = 0.013. In (b) and (d), 
dashed lines represent the temporal profiles of the slow variable. Dashed-dotted 
and dotted lines are the solutions of the reduced equations 1)23(1 - (12411 . Other 
parameter values for FHN are a — 0.1 and b = 0.5. 
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Figure 9: Pulses for FHN with a = 0.1,6 = 0.5 and increasing A. When D = 
0.03: (a) A = 0.01, (b) A = 0.03, (c) A = 0.065. When D = 0.7: (d) A = 0.05, 
(e) A = 0.1, (f) A = 0.237. 
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Figure 10: Pulses for FHN near X C (D) with a = 0.1,6 = 0.5 and increasing 
D. Temporal profiles when (a) D = 1, (c) D = 10, (e) D = 100. The small 
oscillation in the tail indicates that (0, 0) has just become a spiral point. Spatial 
profiles when (b) D = 1, (d) D = 10, (f) D = 100. 



maximum value 2. The width of the pulses does not vanish cither. Moreover, it 
is quite large near the critical coupling D c . Close to X C (D), the peak contains 
3 points if D = 0.009. This number increases to 8 points when D — 0.004. 
For D = 0.003, we reach 30 points. The number of points in the peak of a 
propagating pulse near the critical time scale separation X C (D) seems to increase 
as D approaches D c . Then, why do pulses fail to propagate? Near D c the 
leading edge advances in the slow 'saltatory' manner typical of wave fronts at 
the depinning transition 8 . Then, the separation of time scales described in 
earlier Sections becomes blurred. Both variables, v n and w n change as the 
leading front moves. For D fixed, and as A increases, the leading edge becomes 
pinned and ceases to propagate at a critical value of A. Thus, in the highly 
discrete region, waves fail to propagate at X C (D) due to the pinning of the 
leading edge. 

For intermediate values of D, the peak of propagating pulses near the critical 
time scale separation X C (D) contains one point, see Figure^! (c). The temporal 
profiles of these pulses show a peak with a narrow trapezoidal shape. Thus, their 
spatial profiles change noticeably as the pulse moves, showing the configurations 
with variable height depicted in Figure|Hl(c). The reason for propagation failure 
in this region is the vanishing width of the peaks as A increases. 

For larger D, the mechanism of propagation failure as A grows is related 
to the decreasing height of the pulses, as in the continuum limit. In fact, the 
profiles and speeds of the pulses near X C (D) resemble the profile and speed of the 
pulse solutions of J2J) near the critical time scale separation A c for the continuous 
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problem, as shown by Figures (f) and 1101 As D grows, more and more points 
accumulate in the leading and trailing fronts (see Fig. I10|) and propagation fails 
before the top of the pulse is reduced to one point. 

In this Section, we exploit (and correct) the predictions provided by the 
asymptotic construction of travelling wave trains and pulses carried out in Sec- 
tion^to characterize the critical coupling D c and the critical separation of time 
scales for propagation X C (D). First, we consider travelling pulses and then we 
turn to wave trains. 

4.1 Propagation failure for pulses 

For a fixed D, the asymptotic construction in Section [2 . 31 yields upper and lower 
bounds on the critical value X C (D) above which propagation of pulses fails, as 
a function of D. Moreover, we find an asymptotic prediction of \ C (D) for small 
D, which agrees reasonably with numerical measurements. Prior to this work, 
the only work in this direction is contained in the paper by Chen and Hastings 
who found rather coarse upper bounds on the separation of time scales for 
a particular class of sources that excludes FHN and most realistic sources. 

An upper estimate X C (D) + on X C (D) is found observing that our asymptotic 
construction is consistent when the peak has one or more points: Li(w°, D, X) > 
1. Therefore: 



On the other hand, numerical solutions show that pulses fail to propagate 
when L\(w , D, X) is smaller than a certain number of points k(D). Once k{D) 
is determined: 



yields a lower bound X C (D)~ for X C (D). Without a theory to predict k(D), our 
choice k(D) = 5 for D € (0.01, 1) in Figured (a) is empirical. 

Figure |7fa) compares the numerically calculated X C (D) with the upper and 
lower bounds A c (-D) ± for the FitzHugh-Nagumo model. The numerical approx- 
imation of X C (D) is obtained as follows. We solve the initial value problem 
<|14|) - (|15[l starting from the equilibrium state (0,0) and excite the left boundary 
for a large enough time T: v (t) = 2%(0, T), Wo(t) = 0, t > 0. If a pulse is 
successfully generated, we increase A slightly. The process stops when we reach 
a value of A for which pulses cease to be formed. 

As Figure E£a) shows, bounds i|20|) - (|21|) are not sharp. There are several 
reasons for this. First, the minimum number of points in the peak before prop- 
agation failure changes with D. Second, as A grows the height of the propagating 
pulses diminishes and our asymptotic construction loses precision. Third, near 
the critical coupling D c , the reduced bistable equation does not describe 
accurately the speed of the leading front. The relevance of these three factors 




(20) 




(21) 
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varies with D. We will improve our predictions of \ C (D) by a more precise 
study of the three regions mentioned before: highly discrete, continuous and 
intermediate. 



4.1.1 Highly discrete limit 

In the highly discrete limit, fl < 1, ro° ~ w c -(D) and we are near the depinning 
transition of wave fronts characterized by a saltatory motion of active points, 
see Section 12.11 During the long time intervals between abrupt jumps of the 
saltatory motion, w n can change appreciably and the reduced equation for the 
motion of the leading edge of pulses: 

-^ = D(v n+1 - 2v n + v n -x) + f(v n ,w°) (22) 

is no longer accurate. Instead of (|12fl with w = u>°, the evolution of the active 
point is governed now by: 

dv 

— =D{z 3 {w)-2v + z 1 (w)) + f(v,w), (23) 
dw . , . 

— =Xg(v,w). (24) 

Figures HO (b) and (d) compare solutions of (|23fl - l|24|) . represented by dashed 
lines, with the temporal profile of the pulses. The leading edges are accurately 
described. 

As we increase A, we find a critical value Xq(D) at which trajectories of (|23[1 - 
<|24|) jumping from a neighborhood of (u°, w°) = (zi(w°), w°) to a neighborhood 
of (z 3 (w°),w°) cease to exist. We obtain a first impression about A (D) by 
examining the phase plane in Figure ITT1 System (|23l) - l|24[l has an equilibrium 
point. For A small, the equilibrium point is a node and trajectories starting 
at (zi(w°),w°) reach a neighborhood of (z 3 (w°),w°), as shown by the dashed 
lines in Figure ^] (a). As A grows, it becomes a spiral point. Trajectories 
originating at (zi(u>°), w°) wrap around the equilibrium point, without reach- 
ing (z 3 (w°),w°), see Figure ITT1 (b). If we compute the value of A at which 
the equilibrium changes type, we find an approximate value for \ (D). This 
approximation happens to be quite poor: it is almost constant for D small. 

A better prediction of Xq(D) is found by calculating two characteristic times. 
Suppose w does not depart appreciably from w° « w c - (D). The depinning anal- 
ysis explained in [5] shows that the quiescent period between successive jumps in 
saltatory motion is t v ~ . = During this quiescent period, v ~ sq, 

^/q/3(-u)°-uj c _(_D)) 

so being a value between zi(w°) and z 2 (w°) (see Section 2.1). Suppose on the 
contrary that A is so large that w increases from w° to w c -(D) while v ~ sq. 
This occurs in a time t w = J w c-( D ^ dw/\g(so,w). If t v <C t w , saltatory motion 
of pulses can proceed. In the opposite limit, t v 3> t w , w evolves towards the 
equilibrium value of J23J-IO an d the pulse stops. The crossover value at which 
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Figure 11: Phase plane for the reduced system P3 |l -l|2l jl : (a) A = 0.01, (b) 
A = 0.1. Parameter values for FHN are D = 0.01, a = 0.1, b = 0.5. 



t v ~ t w yields the estimated critical value Xo(D): 

^ ^WE^m r- {D) (25) 

7T Jw° 9{so,w) 

For FHN, w° = and C^ [D) -H- = -Mn(l - b^^) ~ 22s=M. Thus. 

M D ) ^ Wc ~ s ( D) V^i r^- (£>)) . Figured (b) compares A C (L>) with the numer- 
ically measured value of Ao(-D), our asymptotic prediction l|25|l and the upper 
bound X + (D). Now, this upper bound is computed using the approximation of 
the speed provided by instead of the speed of the wave front solution 

of |[BJ|. Notice the good agreement between X C (D) and Aq(-D) for small D. 



4.1.2 Intermediate region and Continuous limit 

In the intermediate region, the peak of the pulses contains one point near X C (D) 
and the upper bound X + (D) provides a good estimate of the critical threshold. 
Formula l|20|l is more accurate when we replace the speed of the wave front 
solutions of the reduced bistable equation (used in Figure [7| (a)) by the speed 
found solving numerically l|23 |) -i|24 p (used in Figure (c)). 

For larger D, the upper and lower bounds X ± (D) yield poor approximations 
due to the decreasing height of pulses. As Figure (a) shows, X C (D) increases 
slowly to the threshold A c for the continuous problem. Notice that equations 
CJ-© with D{h) = -^2 are the method of lines for the numerical approximation 
of the solutions of the continuous system J3J as h — * 0. For large D, the 
continuum limit yields good approximations to thresholds and pulse profiles. It 
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is well known that has two families of pulses: A fast (stable) family and a 
slow (unstable) family. For discrete models with piecewise linear sources, fast 
and slow families of pulses were constructed by Tonnelier . At a critical A 
both families collapse and disappear. This suggests that a similar bifurcation 
might take place at X C (D) for general discrete models. 

4.2 Propagation failure for wave trains 

The asymptotic construction of wave trains in Section 12.21 is restricted to the 
range w £ (w , w c -(D)). For D fixed, wave trains die out asin w c _(D)~ due 
to the pinning of the leading wave front. The speed of the wave train can be 
approximated by the formula c ~ y/ ' a(D)[3(D){w — w c -{D))/ir, obtained in |S] 
by the depinning analysis of wave fronts in discrete bistable equations. When 
w is near the thresholds w c ±(D), the profiles develop steps as shown by Figures 
E]and El Propagation is 'saltatory'. As D decreases, the range (w°,w c -(D)) 
for which wave trains exist becomes smaller, and it disappears at the critical 
coupling D c such that w c -(D c ) = w°. 

Our asymptotic construction also yields estimates on the critical values of A 
above which propagation of trains fails. For each w € (w°, w c ~(D)), there is a 
threshold X c (w, D) below which wave trains with speed c + (w, D) are found. An 
upper estimate X c (w, D) + on X c (w, D) is obtained observing that our asymptotic 
construction is consistent when the peak has one or more points: Li(w, D, A) > 
1. Therefore: 

X < X c (w, D) < X c (w, D)+ = c+{w, D) / (26) 

The largest value of X c (w, D) + is attained for pulses: X C (D) + = X c (w°, D) + . A 
lower bound is found proceeding in a similar way as we did for pulses. 

5 Relaxation oscillations 

When the nullclincs of ©-1211 intersect at an unstable equilibrium (v°,w°) lo- 
cated in the central branch of the cubic, excitability is lost and spatially uniform 
profiles correspond to a stable limit cycle of system (0} for A < A c . Let V(t) 
and W(t) denote this limit cycle. Its period T > 0, increases proportionally to 
A" 1 as A -> 0. Fi guresEta)-(b) illustrate the structure of the limit cycle. 

In this case, periodic wave train solutions of £EJ-© are generated using 
scaled periodic configurations v n (Q) = V(n/c), w n (0) = W{n/c) as initial data, 
for any c > 0. For continuous systems, this fact was pointed out by Maginu 
|2(J|. Figures IT21 (c)-(d) show a periodic wave train for c = 2. We find wave 
trains with spatial period T/c and speed c. In contrast with the excitable case, 
such wave trains cannot be generated from an equilibrium state by periodic 
excitation of a boundary. 

More often, we expect nonuniform profiles v n (t) and w n (t) to display syn- 
chronization phenomena. Figure 1131 illustrates synchronization of a popula- 
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tion of N = 50 oscillators with zero Neumann boundary conditions. Initially, 
half the population is synchronized at a certain phase and the other half at 
a different phase. The trajectories v n {t), w n (t) of individual points behave as 
V(t + <p n ),W(t + 4> n )i with a slowly varying phase <f> n which may become in- 
dependent of n as t — > oo, as shown in Figure IT31 Notice that the diffusive 
coupling is only active in the part of the cycle where v n changes abruptly. 

For a better understanding of the evolution of different initial configurations 
we need an equation for the phases 4> n {t). We will show below that, for small 
couplings D -C I, the oscillator phases <f) n {t) obey the equation: 

= yi ^ V( - 6n + " M ~ 2V ( 9 n) + V(8 n + K-l ~ (27) 

where 6 n = t + </> Tl (r) and (yi(6), y% {9)) is a T-periodic solution of the adjoint 
homogeneous problem: 



yx,t \ ( Mv,w) x 9v (v,w) \(yi\ = (Q 
m J V x 9w (v,w) U 1 o 



(28) 



normalized by ^ J T 2/i(0)V(0) + J/ 2 (0)W(0) = I. For small A, system (J2HJ 

decouples and y\ is computed from yi,t+/ u (V r , W)y\ = with ^ yi{9)V'(6) — 
I. In general, V and are not known explicitly. However, an asymptotic 
approximation is found using the separation of the time scales. The profiles are 
composed of regions where V varies rapidly in the time scale t whereas W is 
almost constant and regions where W varies slowly in the time scale Xt whereas 
V is at equilibrium, f(V, W) = 0. 

Equation (|27|l is obtained as follows. Set r = Dt. Inserting the expansions 
v n = V{t+dp n (T))+Dv£\t,T)+0(D 2 ),w n = W{t+(t> n (T))+Dw^(t,T)+0(D 2 ) 
in CJ-0 we get: 

Dv^\ + (I + D^ T )V = f(V, W) + D[f v (V, W)v$ + f w (V, W)wi 1) ] , , 
+D[V(t + <j> n+1 (T))-2V(t + ^ n (T))+V(t + <j> n ^(r))}+0(D 2 ) V ; 

Dw£\ + (1 + Dcj> n>T )W' = Xg{V, W) + XD[g v (V, W)v { n 1] + g w (V, W)wlP} 



+0(D 2 ) 



(30) 



This yields: 



Vnl ~ [fv(V, W)v { n 1] + f w (V, W)w i n 1) ] = -<f>n,rV'(t + 0„(r)) 

+[V(t + <pn+i{r)) - 2V(t + 4> n {T)) + V(t + 0„-i(r))] (31) 

Wnl ~ %t,(V, W)viP + 9w (V, W)wi 1) } = -4>n,rW'(t + 0„(r)) 

The solvability condition for system (|3I[1 is precisely our equation 1)2 7|). Equa- 
tions for the time evolution of the oscillator phases in continuous problems were 
derived in [231 IT7| . For discrete problems with diffusion in both the fast and 
slow variables see |2"H. 
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Figure 12: For ML: (a) Limit cycle and (b) Periodic trajectory for an uncoupled 
unit; (c) spatial profiles at different times: v n (t) (asterisks), w n (t) (circles); (d) 
temporal profiles of v n (t) (solid line) and w n {t) (dashed line) for n — 20 fixed. 
Parameter values are / — 0.375, D = 0.5, A = 0.01. 




Figure 13: Synchronization in the Morris-Lecar model. We have superimposed 
the time evolution of N — 50 oscillators. Parameter values are I = 0.375, 
D = 0.5 and A = 0.01. 
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6 Conclusions 



We have studied wave propagation and oscillatory behavior in excitable and 
self-oscillatory one dimensional systems with one fast and one slow variables. 
In excitable systems, we have constructed asymptotically a one-parameter fam- 
ily of stable wave trains with increasing periods and speeds. Our numerical 
experiments suggest that the stability of wave trains in discrete systems can be 
characterized analytically in terms of the dispersion law, as it happens for wave 
trains in continuous systems |19| . When the spatial period tends to infinity, we 
obtain the fastest wave of the family: a solitary pulse. As the strength of the 
coupling between nodes decreases, the family of wave trains becomes smaller. 
The disappearance of the solitary pulse marks the onset of propagation failure 
due to spatial discreteness. The critical coupling is the pinning threshold for 
the leading edge of the pulse in a reduced bistable equation. Close to failure, 
the wave profiles develop steps and propagation becomes saltatory. Instead of 
propagating smoothly, the excitation jumps from node to node. Another source 
of propagation failure is the small separation of time scales for excitation and 
recovery, shared by continuous models. We obtain bounds for the critical sepa- 
ration of time scales as a function of the coupling. Technically, we have greatly 
improved the rough bounds contained in , derived under very restrictive as- 
sumptions on the nonlinearities, which exclude FHN and most realistic sources. 
Such bounds are unknown in the continuum case. For small D, we also develop 
an asymptotic theory of failure as A decreases. Our asymptotic predictions show 
a good agreement with numerical solutions of FHN and ML models. 

In self-oscillatory systems, wave trains are generated choosing scaled peri- 
odic configurations as initial data. For other initial configurations, we observe 
synchronization phenomena. In the limit of small coupling, we find an equation 
for the time evolution of the oscillator phases. A detailed analysis of this equa- 
tion which would support the numerical evidence of synchronization remains to 
be done. 
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